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1 Summary 

A new, high-order, conservative, and efficient method for conservation laws on 
unstructured grids is developed. The concept of discontinuous and high-order 
local representations to achieve conservation and high accuracy is utilized in a 
manner similar to the Discontinuous Galerkin (DG)[1] and the Spectral Vol- 
ume (SV)[2] methods, but while these methods are based on the integrated 
forms of the equations, the new method is based on the differential form to 
attain a simpler formulation and higher efficiency. Conventional unstructured 
finite-difference (FD)[3] and finite- volume (FV)[4] methods require data re- 
construction based on the least-squares formulation using neighboring point 
or cell data. Since each unknown employs a different stencil, one must repeat 
the least-squares inversion for every point or cell at each time step, or store 
the inversion coefficients. In a high-order, three-dimensional computation, the 
former would involve unpractically large CPU time, while for the latter the 
memory requirement becomes prohibitive. In addition, the finite-difference 
method does not satisfy the integral conservation in general. By contrast, the 
DG and SV methods employ a local, universal reconstruction of a given order 
of accurac}' in each cell in terms of internally defined conservative unknowns. 
Since the solution is discontinuous across cell boundaries, a Riemann solver 
is necessary to evaluate boundary flux terms and maintain conservation. In 
the DG method, a Galerkin finite- element method is employed to update the 
nodal unknowns within each cell. This requires the inversion of a mass matrix, 
and the use of quadratures of twice the order of accuracy of the reconstruction 
to evaluate the surface integrals and additional volume integrals for non-linear 
flux functions. In the SV method, the integral conservation law is used to up- 
date volume averages over sub cells defined by a geometrically similar partition 
of each grid cell. As the order of accuracy increases, the partitioning for 3D 
requires the introduction of a large number of parameters, whose optimization 
to achieve convergence becomes increasingly more difficult. Also, the number 
of interior facets required to subdivide non-planar faces, and the additional 
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increase in the number of quadrature points for each facet, increases the com- 
putational cost greatly. 

In the spectral difference (SD) method, the conservative unknowns in each 
cell are the number of nodal values required to support a. reconstruction of 
a given order of accuracy. Their locations are chosen so that a quadrature 
approximation for the volume integral exists at least to the same order of 
accuracy. The fluxes are calculated at a different set of nodes, whose num- 
ber will support a reconstruction of one order higher accuracy, since the flux 
derivatives are used to update the conservative unknowns. They are located so 
that quadrature approximations fo r sumacs intocrpic! ov^r r'oll Koii'nde^'i^s 
exist to a required order of accuracy. In addition, the locations of the conser- 
vative nodes and the flux nodes must be such that the integral conservation 
law is satisfied for the cell to the desired order of accuracy. If the nodes are 
distributed in a geometrically similar manner for all cells, the discretizations 
become universal, and can be expressed as the same weighted sums of the 
products of the local metrics and fluxes. These metrics are constants for the 
line, triangle, and tetrahedron elements, and can be computed analytically for 
curved elements. We can also show that the number of flux nodes is far less 
than the number of quadrature points in the SV method. Since all unknowns 
are decoupled, no mass matrix inversion is required. The SD formulation for 
the line element is a generalization of the multidomain spectral method[5]. Its 
tensor products can be used for quadrilateral and hexahedron elements. 


2 The Spectral Difference Method 


The most general form of a conservation law can be written as 


+ V * F = 0. 
dt 


( 1 ) 


where the conservative variable u can be a scalar or a vector, and the general- 
ized flux F can be a vector or tensor. The term V*T represents the divergence 
or curl of F, depending on the physical definition of u. Integrating (1) over 
cell ?', we obtain 


a 

dt 




dS * F = 0, 


( 2 ) 


where V* is the volume of the cell i , and S^i is the area of face l for cell i. (In 
2D, each face is actualty a line.) Here M is the number of faces, which is one 
more than the dimension D of the domain. 

In each cell, the discrete unknowns are the values of u at quadrature points 
for the volume integral over the cell. We denote these points, some of which 
may lie on the cell faces, as i and define Ujj as u(iy ? -). The expansion of u 
in the cell can be written in the cardinal form 



Discontinuous Spectral Difference Method 


An 

u i( r ) = ( 3 ) 

3 = 1 

where Lj^( r) are the cardinal basis functions and N n is the number of basis 
functions required to support a desired degree of precision n of the reconstruc- 
tion. We will use polynomials as an independent basis. The locations of 
then uniquely define the L Ja (r). In order to evaluate the surface integrals in 
(2) efficiently, we discretize F at points most or all of which are located 
at quadrature points for those integrals. The expansion of F in the cell can 
aiso be written in the cardinal form 

A 7 «+i 

F i (r)=^M M ( r)F M> (4) 

k=l 

where Mk.%( r) are now the set of cardinal basis functions defined by r k.i and 
Fk.i = F{vk.i)- We can satisfy (I) at points by evaluating 


V * F(vj :i ) = E VAf*,<(r it< ) 
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( 5 ) 
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In order to evaluate Fk % i , 
(3) as 


Uk.i is required, which can be obtained directfy from 


N n 

= T Lj,i(Tk m i)Uj 


(6) 


3=1 

To reduce the cost of this interpolation, some of points Tk ,» may be chosen to 
coincide with Tj^. If the points Vj % i and i are distributed in a geometrically 
similar manner for all cells and within each cell i the gradient of a function 
is expressed in terms of its area vectors S* and volume Vi, the coefficients in 
(5) and (6) become universal, independent of cell i. There are only a few of 
these coefficients, which can be calculated exactly and stored in advance. For 


points Tk.i in the interior of cell, Fk.i = F(uk.i )■ For those points located 
on the cell faces, since u may be discontinuous, we must replace the flux F by 
a Riemann flux. 

In order to check the integral conservation law for the cell, we must show 


that 


M 


* Ftes) = E s * * E wkFk ’ 

3=1 


(7) 


i-i 


Here Wj are the volume quadrature weights at the points and Wk are 
surface quadrature weights for face Z, and for each face the summation is over 
those points rk,i located on that face. The total contribution from any interior 
point Fk,i to the volume integral in Eq. (7)(righthandside term) must vanish. 

From Eqs. (5) and (6), we also see that the SD formulation is very similar 
to that of the FD method for structured grids. The SD method thus retains the 
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simplicity and computational efficiency of the structured FD method. How- 
ever, the metric terms in the latter are evaluated by numerically differencing 
the grid point coordinates. Since numerical grid generators are mostly only 
second-order accurate, the overall accuracj 7 of the solution can be severely 
degraded if the grid is not sufficiently smooth. In contrast, the metric terms 
in the SD method are computed exactly from the geometry of the grid, no 
matter how it was generated. It thus retains its formal accuracy, even for very 
unsmooth unstructured grids. Furthermore, in contrast to the FD method, 
the integral conservation law is satisfied to the desired accuracy. 


3 Locations of the Unknowns and Flux Points 


The critical part of the SD method is the locations of the u points ly* and F 
points which are determined by symmetry groups associated with the cell 
centroids, vertices, edges, and faces. All but the first contain arbitrary param- 
eters that can be varied to obtain optimum solutions. The number of points 
required for a reconstruction with a specified degree of precision is greater 
than the minimum number of Gaussian quadrature points for that precision. 
One can obtain greater efficiency by locating some u points to coincide with 
F points-: For F points on the vertices in 2D and edges and vertices in 3D, 
more than one Riemann solver is necessar} 7 . For those formulations with ex- 
pensive Riemann solvers, these points should be minimized. Another criterion 
for the placement of u and F points is that the reconstruction matrix is non- 
singular. The final criterion is that integral conservation is satisfied within the 
desired degree of precision. We can show that the number of F points is far 
less than the number of flux quadrature points in the SV method with the 
same accuracy. 

We first show some representative placements of u points (circles) and 
F points (squares) that satisfy integral conservation, for various orders of 
accuracy (which are one more than the degree of precision). Figures la and lb 
show the placements for line element with second and fifth order of accuracj 7 , 
and Figures 2a, 2b, and 2c show' the placements for the triangular element with 
first, second, and third order of accuracy, respectively. Here the u points are 
placed at the Gaussian quadrature points, while the F points on the surface 
are placed at Gauss-Lobatto points, except in Fig. 2a, where they are located 
at the Gaussian points. 


4 Numerical Results 

We first show some one-dimensional convergence studies for the linear con- 
vection equation. Plotted are L\ (solid lines) an (dash lines) error norms 
as functions of grid size for second to fifth order accurate methods. Figure 3a 
(1D-G) shows the plots for the u points at the Gaussian quadrature points, 


Discontinuous Spectral Difference Method 


(a) 2 nd order ' (b) o th order 

Fig. 1. Placement of unknowns and flux points for line element. 



(a) I st order (b) 2r d order (c) 3 rd order 

Fig. 2. Placement of unknowns and flux points for triangular element. 


while Figure 3b (1D-GL) shows the plots for the u points at the Gauss- Lob at to 
quadrature points. Both exhibit the expected orders of accuracy. In Fig. 3c 
(2D) we show analogous plots for a plane wave propagating at 45 degree 
through a square domain for second and third order accurate methods, using 
the point placements in Figs. 2b and 2c. The third example is the scattering 
of a TM wave incident on a perfectly conducting circular cylinder. Figure 
4 shows the unstructured grid consisting of 2024 triangular cells. The wave 
propagates from the left to the right w r ith the wave number equal to 5, based 
on the radius of the cylinder. This gives approximately 6 cells per wavelength. 
Contour plots for E z wTfch the exact solution (solid lines) are shown in Figs. 5a, 
5b, and 5c for first, second, and third order accurate methods, respective^. It 
is seen that the first order solution is very dissipative with this grid resolution, 
while the second and third order solutions show an excellent agreement with 
the exact solution. 
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Fig. 4. Grid for a region exterior to a circular cylinder. 



(a) 1 st order (b) 2 nd order (c) 3 rrf order 

Fig. 5. Contour plots of E z for a plane wave incident on a perfectly conducting 
cylinder, (numerical solutions: color contours, exact solution: solid lines) 





